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1. Introduction 

It is well-known that the dynamics of many physical systems can be derived 
from variational principles. For classical mechanics, for example, the first variation 
of an action functional yields the equations of motion, which are also called the 
Eulcr-Lagrange equations of the functional. By discretizing these equations one can 
obtain numerical methods. This strategy is time-tested and reliable. More recently. 
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2 MICHAEL WESTDICKENBERG AND JON WILKENING 

however, an alternative approach has attracted considerable interest, which consists 
of discretizing the action functional instead. The numerical method is then given by 
the Euler-Lagrange equations of this discrete functional. Methods obtained in this 
way enjoy remarkable stability properties and preserve the "symplectic" structure 



of the problem. We refer the reader to 13 for further information. In the present 



paper, we apply a similar strategy for the porous medium equation and the system 
of isentropic Euler equations in one space dimension, and derive numerical methods 
by discretizing a suitable variational principle. 

The isentropic Eulcr equations model the dynamics of compressible fluids under 
the simplifying assumption that the thermodynamical entropy is constant in space 
and time. These equations form a system of hyperbolic conservation laws for the 
density g and the Eulerian velocity field u. They take the form 

dtg + V -{911)^0, 

dt{gu) + V • (£-u ® u) + VP((?) = 0, ^ ' ' 

where {g, u) : [0, 00) x R'^ — > V are measurable functions with U := [0, 00) x M.'^. 
Since entropy is assumed to be constant, the pressure P{g) depends on the density 
only. It is defined in terms of the internal energy U{g) of the fluid by 

P{g) := U'{g)g - U{g) for ah g^Q. (1.2) 

For the important case of polytropic fluids, we have 

Uig) = -^ and P(g) = Kg"^. (1.3) 

7- 1 

Here 7 > 1 is the adiabatic coefficient, and we will assume the common normaliza- 
tion K := O'^/j with :— {j — l)/2. For isothermal fiows, we have 

Uig)^g\ogg and Pig) ^ g. (1.4) 

We are interested in the Cauchy problem, so we require that 

(e,u)(t = 0, •) = (^,u) 

for suitable initial data [g, u) with finite mass and total energy. 

It is again well-known [2 10 that the system of isentropic Euler equations can 



be considered as the Euler-Lagrange equations of a certain action functional, so we 
might try to derive numerical schemes by following the procedure outlined above. 
But we want to work in the framework of optimal transport theory. Our discussion 
in the present paper will therefore not be based on the variation of an action 



functional, but on the notion that the system of isentropic Euler equations ( 1.1 ) can 
be considered as a steepest descent on a suitable abstract manifold. This notion was 
introduced in [8] and will be explained in detail below. To motivate our approach. 



recall that any smooth solution (p, u) of (1.1 1 also satisfies 
1 



dt\^-g\u\' + U{g)j+V-U-g\u\' + Qig)ju\ =0, (1.5) 

where Q{g) ;= U{g) + P{g). Therefore the functions 

i]{g,u):=-g\u\^ + U{g) and q{g,u) -.^ (-g\u\^ + Q{g)ju 
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form an entropy /entropy flux pair for ( l.ll; see |6i. The entropy ri{g, u) is a convex 



function of g and gu. As a consequence of (1.51, the total energy 



E[g,u]it):= f Ug\u\^ + U{g)){t,x)dx (1.6) 

is conserved in time as long as the solution remains smooth. 

On the other hand, it is known that solutions of ( |1.1[ ) are generally not smooth: 
No matter how regular the initial data is, jump discontinuities can form in finite 
time. These jumps occur along d-dimensional submanifolds in space-time and are 
called shocks. Across shocks, total energy is dissipated (that is, transformed into 
forms of energy such as heat that are not accounted for in an isentropic model) , so 



(1.5) cannot hold anymore. It is therefore natural to consider weak solutions (g, u) 



of ( 1.1 1 such that E\g^ u] is a non-increasing function of time. 

As will be explained below, we consider solutions of the isentropic Euler equations 
as curves on the space of probability measures, using optimal transport theory. We 
are interested in curves along which the total energy decreases as fast as possible, 
under conservation of mass and momentum. This steepest descent interpretation 
has already proved very fruitful for certain degenerate parabolic equations, such as 
the porous medium equation; see below. Numerical schemes based on this varia- 



tional principle have been derived in 11 . As the authors point out, these schemes 



have remarkable properties: oscillations can be reduced because the optimizations 
take place in weak topologies; the approximations can be constructed using dis- 
continuous functions because no derivatives are needed; and overall, the schemes 
are very stable. In this paper, we revisit the discretization of the porous medium 
equation, and also introduce variational approximations for the isentropic Euler 
equations. We demonstrate that our schemes capture very well the nonlinear fea- 
tures of the flows. The simplest versions of our schemes are first order accurate. 

Designing higher order methods in the optimal transport setting is an interesting 
challenge that has not previously been addressed. We present schemes that achieve 
second order convergence (away from shocks) and propose a general strategy for 
constructing higher order methods in a backward differentiation formula multi-step 
(BDF) or diagonally implicit Runge-Kutta framework. Again we observe that our 
schemes, due to their variational nature, are remarkably stable. 

To put our work and that in l8| in perspective, let us first recall recent research by 
various authors on a steepest descent interpretation for certain degenerate parabolic 
equations. It was shown by Otto fl5' that the porous medium equation 

Sfg- AF(£.) =0 in[0,oo)xM'^ (1.7) 

is a gradient flow in the sense enumerated below. This result has been generalized 
considerably, and we refer the reader to [1 17 for further information. 



(1) We denote by ^{M.'^) the space of all £''-measurable, non-negative functions 
with unit integral and finite second moments, where £'' is the Lebesgue measure. 
The space ^(M'*) is equipped with the Wasserstein distance, defined by 



W{gi,g2f:=mil // |x2 - xip 7(^X1,^X2) : ^^#7 = ft-C'^ >, (1-8) 




where the infimum is taken over all transport plans 7, which are probability mea- 
sures on the product space K'^ x K'' with the property that the pushforward 7r'#7 
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of 7 under the projection tt* : M'' x M^ — > R'^ onto the ith component equals gi2'^. 
The infinuim m (1.81 represents the minimal quadratic cost required to transport 



the measure £»i£'' to the measure g2^'^, and one can show that it is always attained 
for some optimal transport plan 7. In fact, there is a lower semicontinuous, convex 
function (p: M'^ — > M, uniquely dete rmined gi-a..e. up to a constant, such that 
7 = (id X V(/))#(eii£''); see 0[3)[7)[l7]. This identity implies that 



W(ei,e2)'=/ \W<l>{x)~x\''gi{x)dx and (V(/.)#(ei£'^) = ^2^''. 

JM'' 

We call V0 an optimal transport map. It is invertible gi-a.e. and the Hessian D^(j) 
exists and is positive definite gii-a.e. When d = 1, the Wasserstein distance and the 
transport map can be computed more explicitly; see below. 



(2) We introduce a differentiable structure on 



as follows: For any point 



£•€ 



), the tangent space T^^ 



) is defined as the closure of the space of 



smooth gradient vector fields in the ^^(M'*, £))-norm. This definition is motivated 
by the fact that for any absolutely continuous curve 1 1-^ gt E ^{M.'^) with go = g, 
there exists a unique u S Tg^(E'^) with the property that 



^*^*lt=o 



V-{gu)^0 in ^'(M"^). 
That is, for all test functions (p G S^{R'^) we have 

gt{x)(f){x) dx — I u{x) ■ '^(t){x)g{x)dx, 



(1.9) 



d 
dt 



t=o 



and so a change in density can be accounted for by a flux of mass along the velocity 
field u. This structure makes :^{R'^) a Riemannian manifold, and we define 



T^(I 



U{(^'") 



ge 



,ueT 



„^iR'^)\. 



(3) li u = —g ^VP(gi), then ( |1.9[ ) yields the porous medium equation (1.7| at 
one instant in time. This vector field is the "gradient" of the internal energy 



U[g] 



U{g{x)) dx with P{g) = U'{g)g - U{g) 



(1.10) 



in the sense that u is the uniquely determined element of minimal length in the 
sub differential oiU[g] with respect to the Wasserstein distance. The vector field u 
is indeed tangent to ^{M.'^) because —g^^'^P{g) = —'^U'{g) is a gradient. 

The internal energies for the porous medium equation that arise most frequently 
in physical applications are identical to those of the isentropic Euler equations. 



namely (1.31 or (1.4|. For the latter choice with k = 1, equation (1.7) becomes 



at£»-A£i = in [0, 00) X M'*, 
i.e., the porous medium equation reduces to the heat equation. 



(1.11) 



The interpretation of (1.7) as an abstract gradient flow suggests a natural time 



discretization: Given a time step r > and the value p" S ^(M'^) of the approxi- 
mate solution at time t" := nr, the value at time i"+^ is deflncd as 



„n-f-l 



argmm • 



2t 



W{g'\gf+U[g]: g d 3^{W^) 



(1.12) 
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We assume that p" has finite internal energy U[g'^] < oo. It can then be shown that 
for all n ^ 0, problem (1.12) admits a unique solution with U[g'"'~^^] ^ UIq^]. We 
describe the motivation for the minimization (1.12) in Appendix [A[ 

It is useful to reformulate the porous medium equation ( 1.7 ) in Lagrangian terms. 



We fix a reference density g and consider a time-dependent family of transport maps 
X{t, ■):R'^- — > W^ such that X[t, ■)i^g = g{t, O^'' for all t ^ 0. If 

X{t, x) = X{t, y) implies that dtX{t, x) = dtX{t, y) 

for g-Si.e. x,y ^ M"^, where d^ denotes the forward partial derivative in time, then 
we can define an Eulerian velocity u by the formula d^X =: uo X. For the porous 
medium flow, the velocity is given by Darcy's law u = — V[/'(£'). This vector field 
is related to the first variation of the functional V[X] := U[X^g], which is defined 



for all £i-measurable maps X : 
in Appendix [Xj we then have 

"^ V[X + eiC o X)] - 

e=0 



Let 



5VIX] 

sx 



VU'{g) o X. As explained 



de 



SV[X] 

sx 



for all smooth test functions ^ : 



■{CoX)dg^ I \/U'{g)-Cgdx 
RecaU that X#g = gZ"^. 



The minimization (1.12) can now be rewritten in terms of transport maps as 



X 



n+l 



argmm 
X 



m|^£|X-X"|^d^ + Vm|, 



(1.13) 



where X" is chosen such that X^#g = g"£'^. We obtain X"+'^ = r"+i o AT" for all 
n ^ 0, where r"+^ is the optimal transport map pushing g^£'^ forward to gi"+^£''. 
The discrete analogue of Darcy's law is then u"+^ = —VU'{g"^^), where 



r"+i - id 



u 



(j.«+i)-i ^ 



a:"+i - X'' 



(X"+i)- 



(1.14) 



We refer the reader to Appendix \K\ for details. In particular, we have that 



= -vc/'(g"+i)oA:"+i = 



sx 



X=X"+i 



which shows that the minimization problem (1.13) is equivalent to one step of the 



_ 5V[X] 



backward Euler method for the dynamical system dtX — jj^ 

Notice that if for each timestep we reset the reference measure g to be equal to 
£("£'*, then we may choose X" = id and obtain Ar"+^ ~ r"+^, so the minimization 



(1.13) reduces to (1.12) if the latter is rewritten as a minimization problem for the 
optimal transport map pushing g^£,'^ forward to £>""'"^£'^; see again Appendix [A| 
This is what we do in our first order schemes detailed below. In designing more 
complicated multi-step or Runge-Kutta schemes, however, it is useful to describe all 
the intermediate steps or stages using one single reference density, which requires 



the generality of the second half of (1.14). We also emphasize that in contrast to the 
space of probability measures ^(M'*), the space of transport maps X : W^ — > K'^ 
is flat, with a linear structure and a translation invariant metric. 

A similar steepest descent interpretation as for the porous medium equation can 



be given for the system of isentropic Euler equations ( 1.1 ). The philosophy here is 



that among all possible weak solutions of ( 1.1 ) we try to pick the one that dissipates 
the total energy (1.6) as fast as possible, as advocated by Dafermos [6]. Note that 
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unlike in the porous medium case, for (1.1 1 the energy dissipation wih be singular 



in the sense that it only occurs when the solution is discontinuous. 

At the heart of the variational time discretization introduced in [8] lies an opti- 
mization problem similar to ( 1.12 ), but with a different homogeneity in the timestep; 



see (1.171 below. In fact, the relaxation of internal energy is a second order effect 



only. What dominates the flow is the transport of mass along the integral curves 
of the velocity field. Let us recall the time discretization introduced in r8j. 

Definition 1.1 (Variational Time Discretization). 

(0) Let initial data {g°, u°) G T^(M'*) and S > he given. Put f := 0. 
Assume that g" € ^(M'') and u" e ^^{R'^, g^) are given, where n ^ 0. Then 

(1) Choose r £ [5/2, 5] in such a way that the push-forward measure 

(id -H ru")#((?"£'') =: e"£'' 
is absolutely continuous with respect to the Lebesgue measure. Find 

(id-t-Tu")#(g"£'^) = ^"£'* 
e^2(M'',e") such that 



u 



(1.15) 



That is, find the velocity field u" with minimal ^^(M'', p")-norm such that the 
map id-|-ru" pushes the measure g^Sf^ forward to g^Sf^. Such a vector field always 
exists and is uniquely determined. It is given by the gradient of a semi-convex 
function, and is therefore a tangent vector. Moreover, we have the estimate 



u"r£»"da;< 



u|^e"da; 



(1.16) 



for any velocity u satisfying the constraint that (id -I- ru)#g" — (p . This shows 
that by replacing u" by the optimal transport velocity u", we decrease the kinetic 
energy as much as possible. The map id + tu" is invertible g^-a.e. 



(2) Update the density by computing the minimizer 



„n+l 



argmm • 



4r2 



W{g'\gf+U[g]: gE ^{R'^) 



(1.17) 



One can show that the density p""*"^ is uniquely determined and that 

2r2 

IT 

2 



Here U[g\ is the internal energy functional from (1.10 1, which is a part of (1.61 

miquely determine 

V[/'(e"+i))#e"+i = ^" 



id 

2t2 



VC/'(e"+i) 



^n+l 



dx 



W((?"+\e")2. 



The density e"+i is thus regular in the sense that ^-VU'{g'"+^) G ^^( 
In fact, the latter vector field is the gradient of a semi-convex function. 

(3) Update the velocity by defining 

2r2 



,9' 



,"+1^ 



,ri+l 



(id 



id + ^-VC/'(g"+i) - rVC/'(e"+i) 



(1.18) 
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(4) Let t"~^^ := t" + r, increase n by one, and continue with Step (1). 

We refer the reader to fsl for more details. It is shown there that 

E[g"+\u''+^] sC E[g'^, u"] for aU n > 0. 

One can then define the function 

{g,u){t,-) :=(e",u") foralHe [e,e+^) and n > 0, 

which is piecewise constant in time and approximates a weak, energy-dissipating 
solution of the isentropic Euler equation (1.1). A formal argument is given in [s] 



for the convergence of this approximation when the timestep (5 — > 0. 

We give a derivation of the time discretization of Definition |1.1| in Appendix [B] 
From the point of view of numerical analysis, it is interesting to generalize this 
scheme to a family of first order implicit methods that also includes the backward 
Euler method as a special case. To this end, let us temporarily assume the solution 
is smooth and write the isentropic Euler equations in a Lagrangian frame as 

1(f) = (rm) ""■ fW:=-^.-W/'(„,oX, ,1.19) 

where X^g =: gxSf^ and V := u o AT is the pull-back of the Eulcrian velocity. 
Again we set V[X] := U[X^g], as in the porous medium case discussed above. Let 
now a € (0, 1] be a fixed parameter and consider the implicit method 

y^+ij = {^V'^)^'^\ f[A"+i] )' ^^'"^^^ 

Using the first equation to eliminate T/"+^ in the second, we obtain 

1 /^«+l _ (^« ^ ^yn^\ ^ f rj^n+ll^ n 21) 

yn+l ^ynj^ rf [A"+i]. (1.22) 



Equation (1.21) will be satisfied by the solution of the minimization problem 



A"+i = argmin i -^ / |X - (A" + Ty")|2 dg + V[A] I. (1.23) 



Once this is solved, we substitute the left-hand side of (1.21) for f[A"+ ] in (1.22) 
to obtain V""^^. This is the Lagrangian reformulation of the energy minimization in 
Step (2) of Definition IlT] As explained there, we combine this step with a velocity 
projection, so we replace in each step the velocity V" by an optimal transport 
velocity, which in the Lagrangian framework takes the form V" = u" o A" with 



u" given by (1.15). This modification dissipates the maximum amount of kinetic 
energy without changing the convected distribution of mass. 

The three choices of a that seem most natural are a = 1, which corresponds to 
the backward Euler method; a = |, which corresponds to the original time dis- 
cretization proposed in [8 (see Appendix B for further information); and a = jj 
which causes the acceleration A"+^ = f [A~^] to agree with its Taylor expansion, 
A"+i « 2r-2(A"+i - (A" + tF")). All three of these variants yield first or- 
der methods that appear to capture shocks and rarefaction waves correctly in our 
numerical experiments. We present a second order version in Section [2?2l 
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2. Variational Particle Scheme 

We now introduce a fully discrete version of the variational time discretization 
discussed in the previous section for the one-dimensional isentropic Euler equations 

dtg + dx{gu) = 

dt{gu)+dAgu^ + P{g))^0 

as well as for the one-dimensional porous medium equation 

dtg~AP{g)=0 in[0,oo)xM. (2.2) 



2 , r./ NN ^ ™ [O.OO) Xl^' (2-1) 



The pressure P is related to the internal energy by equation (1.2), with energy 



density U given by either (1.3) in the case of polytropic gases (we put k = 1 for 



the porous medium equation), or by (1.4) in the case of isothermal gases 



2.1. First Order Scheme. In the simplest version of the algorithm, we assume 
the fluid to be composed of N particles of equal mass m :— 1/N, where A^ > 1. 
At a given time i", the particles are located at positions x" = {a;", . . . ,a;^} with 
velocities u" = {m", . . . ,u^}. In the following, we will always assume that the x" 
are ordered and strictly increasing. The total energy is then given by 

N N-1 , s 

i;(x»,u"):=^-mKp+5:C/ 3^^ (xr+i-xD. 

In accordance with thermodynamics, we think of the density as the inverse of the 
specific volume, which we define as the distance between two neighboring particles 
divided by the particle mass. This formula is used in the internal energy. 

Definition 2.1 (Variational Particle Scheme, Version 1). 

(0) Let initial positions x'', velocities u°, and r > be given. Put t^ :— 0. 
Assume that positions x" and velocities u" are given for n ^ 0. Then 

(1) Compute intermediate positions defined by y^ := x" -I- rw" for all 1 ^ i ^ A^. 
Find the permutation a of the set {1, . . . , N} with the property that the x" :— ya(i) 
are nondecreasing in i. We assume that if there exist indices i < j with Xi — Xj, 
then a{i) < (j{j)- Then define wf := (i" — xf)/T for all i. 

(2) Define new positions x"+^ as the minimizcr of the functional 

N „ N-1 , s 

Z — 1 l — \ ^ ^ 

over all families of positions z e M^ such that Zj+i > z^ for all i. 

(3) Compute the new velocities u"+^ := u" -f ^(a;"+^ — x") for aU i. 

(4) Increase n by one and continue with Step (1). 

Remark 2.2. In Step (1) we perform a free transport of all particles in the direction 
of their given velocities. We obtain a new distribution of mass that is characterized 
by the new particle positions y^. Note that we defined the density in terms of the 
distance between neighboring particles. Since the original ordering can be destroyed 
during the free transport step, we rearrange the new positions in nondecreasing 
order. On the discrete level, it is not necessary to make sure that positions do 
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not coincide because we never compute the internal energy of this intermediate 
configuration. Therefore we can use the same timestep r for aU updates. 

To understand the definition of the new velocities in Step (1), notice that the 



Wasserstein distance as defined in (1.8) can easily be generalized to pairs of proba- 
bility measures /ii and /i2 with finite second moment; see [l] |17| for example. If the 
two measures are given as convex combinations of Dirac measures: 



1 ^ 1 ^ 

where Xi G K'' and yi £ M'' for all i, then the squared Wasserstein distance between 
the measures fii and ^2 coincides with the minimum of the functional 

1 ^ 
W{<7)':=-J2\'^,-y^^^,^\' 

among all permutations a of the index set {!,... ,./V}. This amounts to a Linear 
Assignment Problem, which in the one-dimensional case reduces to sorting the 
positions yi. This can be done with worst-case complexity ©(iVlogTV). In our 
case, the yi are actually largely ordered already, so the sorting is quite inexpensive. 
Step (1) therefore amounts to computing the Wasserstein distance between the old 
and the intermediate positions of particles, and then finding the optimal velocity 
that achieves this transport: We have that x" -I- tm" = x" for all i. 



The functional _F in (2.3) is the analogue of the one in (1.17). Since we minimize 



over positions Zi that are increasing, the first term of (2.3) can again be interpreted 
as a Wasserstein distance squared between convex combinations of Dirac measures 
located at positions Zi and xf . Note that the choice Zi :— x" for all i is admissible 



in (2.3) and results in a finite value of the functional. Therefore the new positions 
x"'^^ must have finite internal energy and thus be strictly increasing. 

Remark 2.3. With only minor modifications we obtain a particle scheme for the 



porous medium equation (2.2), which contains the heat equation as a special case. 



Since now the only unknown is the density g, Steps (1) and (3) in Definition 2.1 
become irrelevant. In Step (2), we substitute the factor l/(2r) for 3/(4r2) to match 



(1.12 ), and the intermediate positions Xi must be replaced by Xi. Of course, we can 
also implement the backward Euler method for the isentropic and isothermal Euler 
equations by changing 3/{4t'^) and 3/(2r) to 1/{2t'^) and 1/r, respectively. 

2.2. Second Order Scheme. In this section, we derive a variant of the scheme of 



Definition 2.1 that is formally second order in space and time. Before doing so, let 
us collect a few facts about the Wasserstein distance in one space dimension. For 
any probability measure /i in M we define the distribution function 

F^(i) ::=^((-oo,t)) foralHeM. 

The function i^^ is nondecreasing, but not necessarily strictly, and can be discon- 
tinuous if the measure ^ contains atoms (Dirac measures) . We define 

F-\s) := sup heR: Ff,{r) < s\ for all s e [0, 1], 

which is the generalized inverse of F^. For any pair of probability measures {fi, v) 
on M with finite second moment, the Wasserstein distance between ji and v can 
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then be computed in terms of the inverse distribution functions as 

W(m,^)':= / \F-\s)~F-\s)\'ds. (2.4) 

If /i does not contain any atoms, then i^^ is continuous, and the composition F~^oF^ 
is an optimal transport map pushing /i forward to v] see Theorem 6.0.2 in [11 . Note 
that the inverse function F~^ pushes the measure l[o,i]£^ forward to ^. 

Assume now that the measure /j, is absolutely continuous with respect to the 
Lebesgue measure, with a piecewise constant density: Let numbers xq < . . . < xn 
and TTii ^ be given such that X]i=i rrii = 1. Let /i = g£^ with 

'O ifx<xoorx>a;^, 

^(^) - ^ —^ if X e [x,_i, x,) andl^is^N. ^2-^) 

i Xi — i 

Then the inverse distribution function of /i is a piecewise linear function that can 
be computed explicitly. Defining Sk '■— X]i=i "^j ^'^^ ^ k ^ N, we have 

Pn ('V == ^k-1 1- Xk tor all s G [sfc_i, SfcJ 

Sk- Sk-l Sk - Sk-l 

and 1 ^ k ^ N. If the second measure v is piecewise constant as well, it is possible 



to compute the Wasserstein distance between /i and v using formula (2.4 1. This 
is particularly simple if both fj, and v have the same number of intervals, and if 
the same mass is assigned to each interval, because then the inverse distribution 
functions i^^^ and F^^^ are piecewise linear on the same intervals. 
Assume now that v = g2^ with 

ii X < xq or X ^ xn, 

if x e [xi^i, Xi) and 1 ^ i ^ N, 

for suitable numbers xq < . . . < xn and rfii ^ with X]i=o ™' ^ ^' ^^ want to 
project ly onto the space of measures of the form ^ = gZ'^ ^ where the density g is 



given by (2.5) and the masses rrii are fixed. More precisely, we want to choose the 



positions x := {xq, . . . ,xm) in (2.5) so that the Wasserstein distance 
W{n,vf^ I \F;\s)-F-\s)\^ds 

-'[0,1] 

is minimal. This amounts to a quadratic minimization problem that can be solved 
easily. Let {</3fe}^o ^® ^^^ standard finite element hat functions with vertices Sk 
as defined above. Then the minimizer is given by x = A^^b, where 

Aki-= / ipk{s)(pi{s)ds and bk ■= / (pfc(s)F~\s) ds 

"'[0,1] "'[0,1] 

for all ^ fc, / ^ N. The calculation of Aki and bk is straightforward. Notice that 
the function F^^^ is piecewise linear on the intervals [sk-i,Sk] with Sk ■= X)i"=i ^i: 
and the Sk are different from Sk (otherwise there is nothing to do). The computation 
of bk therefore requires a partition of [0, 1] that uses both sets of positions. Even if 
the Xi are nondecreasing, the components Xi of the minimizer may not be. 

The projection outlined in the previous paragraph works also in the degenerate 
situation where i^-i = Xi for one or more i. In this case v has a Dirac measure 
at position Xi and rhi denotes the mass of that measure. The inverse distribution 
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function Fj^^ is then no longer continuous, but it is still straightforward to compute 
the numbers 6j. The minimizer x can be obtained as above. 

We can now define our second order method for the isentropic Euler equations. 
Let us first discuss the accuracy in space. Instead of working with point masses, we 
now approximate the density by a function that is piecewise constant on intervals of 
variable length. Specifically, we fix numbers mi ^ with X]i=i "^» ~ 1- ^^ ^'^y time 
i", the density g" is then determined by a vector x" = {xq, x", . . . , a;^) G M^+-'^ of 
positions with x"_2 < a;" for all 1 ^ i ^ A^, through the formula 

if a; < Xq or a; ^ a;^, 

s"(^)--^< "^^ \ixeK_„xf)andl^z^N. (2.6) 

We approximate the velocity at time i" by a piecewise linear function on the inter- 



vals [xf_i,xf). It is determined by u" = (uq , u", . . . , u^) e M^+^ as 



if a; < a;Q or X ^ a;^. 



u"(x) := < x'2 — x 



X- 



■<-i + ^< if ^ e [x«_i,xf) and 1 s^ * s^ iV. 



For any ti ^ consider now the push- forward measure v :~ (id + tm")#(p"£^). 
Let Hi :— x" -I- ru" for all ^ i ^ iV, and let cr be a permutation of {0,1,..., N} 
such that the positions Xi :— ya{i) are nondecreasing in i. We assume that if 
there are indices i < j with Xi ~ Xj, then a(i) < (j{j). The measure v is piecewise 
constant on intervals [xi_i, Xi) for which Xi-_i < x^, and we denote the mass carried 
by this interval by rhi and the density by Qi := rhi/{xi — Xi_i). If Xi_i — x^, then 
the measure v has a Dirac measure located at position Xi and rhi denotes the mass 
of this measure. By construction, we have rhi ^ and X]j=i irii — ^■ 

To compute the numbers rhi we first initialize them to zero. We then run through 
the original intervals [xi_i,Xi) and assign their masses to the new intervals they 
are mapped onto, in proportion to the lengths of these new intervals: 



for 


i = lto N 










k — niin{o-^^(i — 


1) 


cr' 


-Hm 




I = max{cr~^(i — 


1) 


a' 


-'m 




for j = k + 1 to 


I 








if x^ = x^ 










mj = 


rhj 


+ 


rrii 




l-k 



(point mass case) 



else 



rhj — riij + rrii^ ^ — (distribute according to length). 

As we explained above, the computation of the Wasserstein distance in one space 
dimension is very simple if the pair of measures is induced by two piecewise constant 
densities that have the same number N of intervals, and that assign the same mass 
rrii to the ith interval. We therefore project the push-forward z/ = (id-f ru")#(g"£^) 
onto the space of measures of this form, before proceeding with the variational time 
discretization (VTD) algorithm of Definition] 1.1 [ The minimization in ( 1.12 1 is only 



performed over densities g of the same form; see below for more details. Otherwise 
this version of the fully discrete algorithm is identical to the VTD algorithm. 
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Next we show how to modify the algorithm to achieve second order accuracy 
in time. Consider again the Lagrangian formulation (1.19) of the system of isen- 



tropic Euler equations. The key feature of the scheme (1.20) that allows X"+^ to 



be determined by solving a minimization problem of the form (1.231 is that the 
functional f is only evaluated at X"+^ in (1.20). It is also desirable to choose an 
i-stable method to prevent oscillations from growing near shocks. The second or- 
der backward differentiation formula method 9i satisfies both of these conditions. 
Therefore it is a natural candidate to try to adapt to our variational particle scheme 



framework. Applied to problem (1.19), BDF2 takes the form 



yn+l 



X' 



Xn-1 

yn-l 



( V" 



Solving for X :— X"+^ and F"+^ in turn, we obtain 



(2.7) 



(x - (X" + rT/")) -\(x- (X"-i +2t(|V" + IV"^^))] 



f\X] 



T/"+i = 2V'' 



(|F" + |V"""i) + ^f[X'^ 



(2.8) 



The first equation characterizes the solution of the minimization problem (2.9) 



below, which is similar to (1.23). As before, we handle the formation of shocks by 



dissipating the maximum amount of energy possible without changing the convected 
distribution of mass from i"~^ to i" to i"+^. Notice that the convex combination 
|T^" + ^V^^^ (rather than V"'~^) is used to transport mass from time i"~^ to time 
i"+^. As V and u are related via ]/ = u o X, we will use u to denote velocities 



below. Here is our second order method for the isentropic Euler equations (1.1): 



Definition 2.4 (Variational Particle Scheme, Version 2). 
Fix once and for all masses rrii ^ with X]i=i Tfii = ^■ 

(0) Let initial positions x*^, velocities u°, and r > be given. Put t^ :— 0. 
Assume that positions x" and velocities u" are given for n ^ 0. Then 

(la) Compute the intermediate positions yi := x" + ru" for all ^ i ^ A^. 
Sort them and redistribute the masses to obtain i" and rhi. Project back onto the 
space of densities that are constant on N intervals and assign mass nii to the ith 



interval. Let x = [x'q, . 



'-N. 



denote the positions of the minimizer obtained from 



the quadratic minimization problem induced by this projection. See the beginning 
of this section for more details. Then define the velocities u[ := {x'^ — x^)/t. 



(lb) Compute the intermediate positions yi 



2Ta 



imr' 



) for 



all ^ i ^ A^. Sort them and redistribute the masses to obtain x" and rhi as 
before. Project back onto the space of densities that are constant on A^ intervals 
and assign mass rrii to the ith interval. Let x" — (xp , . . . , x'j^) denote the positions 
of the minimizer obtained from the quadratic minimization problem induced by 
this projection. Then de fine the velocities u" := (.x" — x"~"'^)/(2r), which play the 



role of ^V' 



-yn 



in (2 



Note: Skip Step (lb) on the first iteration 
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(2) Minimize the functional 

F{z) := :^^|lz-x'|i:„ - 



2t2' 



8r2' 



Z X 



N 



(z, - z,_i) (2.9) 



over all z — {zq, . 



. , zn) such that Zi_i < Zi for all 1 ^ z ^ A^. We define 

2 



JV 

E 

j=0 



Zi<y5i(s) 



ds, 



where (/3i is the piecewise linear function satisfying ipi{sk) = Sik for all ^ fc ^ iV, 



with Sk ;= J2i=i ' 



In fact, we have ||z 



Hj 



[0,1] 



: Z"^ylz with the matrix A defined by 

ip^{s)ipj{s)ds 



for all ^ i,j ^ N. We denote by x"+^ = (xq'*'^, ■ ■ ■ ,x^^) the minimizer of this 



convex optimization problem. Notice that the Euler-Lagrange equations of (2.9) 



are of the form (2.8) in the second order backwards differentiation scheme. 



Note: On the first iteration, minimize instead the functional 



F{z) := 



4r2 



N 

E 



u 



Zi-1 



(Zi - Zi_l). 



(3) Compute the new velocities 



,n+l 



2«^ - < + ^ (<+i - x^) - 1 (x^ 1 - x^') . 



for all ^ z ^ iV. Notice that this formula is of the form (2.8), where we recall 
that u'l represents |u" + |m"~^, not u"~"^. 
Note: On the first iteration, define instead 

<+ ^ := -^ + |(-r^ - ^ 

(4) Increase n by one and continue with Step (1). 

Remark 2.5. While the £" are nondecreasing by construction, the projection can 
cause the components of the minimizer to be slightly out of order. This is acceptable 



as the minimization (1.17) over densities g of the form (2.61 gives the same result 



if we replace the push-forward measure v := (id + Tu")^{g"£}) by its projection 



onto the space of densities of the form (2.6). Thus, in spite of the projection, we 



still minimize the distance from i' to compute g""*"^. The values a;" will again be 
strictly increasing since otherwise the internal energy would be infinite. 

Remark 2.6. As before, the same procedure works for the porous medium equation 
if we eliminate Steps (1) and (3) and replace the functional F in step (2) by 

N 



F{z) := -||z-x 

T 



rill 2 



h 



,n-l|i2 



E^ 



[Zi 



v). 



On the first iteration, one should use 



F(z) 



^1 



,n||2 



N 

E 

1=1 



u 



Vfli 



[Zi 



Zt-l) 
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Remark 2.7. We have also constructed second (and higher) order variational parti- 
cle schemes in the diagonally implicit Runge-Kutta framework p\. The underlying 
Runge-Kutta method must be "stiffly accurate", i.e., the numerical solution at t""^^ 
must agree with the last internal stage of the Runge-Kutta procedure. We illustrate 
the basic idea using the two-stage second order DIRK scheme with Butcher array 

1/4 1/4 
1 2/3 1/3 . (2.10) 







2/3 1/3 


need to solve is 






+ 4 [i[X-+^) 




2t / T/"+i ^ 



yn+l 

f[x"+i: 



(2.11) 



(2.12) 



Equation (2.11) is simply the backward Euler method with steplength t/4. There- 



fore it can be solved as in ( 1.21 )-( 1.23) above with a, r, X"+^ and y"+^ replaced 
by 1, r/4, X"+4 and V""^^, respectively. Next we choose a suitable linear combi- 
nation of ( |2.11 1 and (2.12) to eliminate the term f[X"+i]. We find 

yn+l] ^\yn+i]+^\yn] -^\f[x"+i 



(2.13) 



The scheme coefficients (2.10) were chosen to minimize the magnitude of the coef- 
ficient 8 in (2.131. Since X"+i and V"^^ are already known from solving (2.11), 



this equation is structurally identical to the multi-step scheme (2.7). As before, the 
unknown X :— X"+^ can be found independently from F""*"^. 

24 " 



We have 



(x - (x"+i + f y"+i)) - -(x - (X" -F (f y" + f v^"+j))) 



fW, 



V 



n+1 



8( 



3y"+j 



) - 5(iy" + jV'+i) + f f[X"+i]. 



The first of these is solved by a minimization problem similar to (1.23) and (2.9) 
above. The left-hand side is then substituted for f [X"+^] in the second equation. 
We leave the details of spatial discretization and energy dissipation (through a 
choice of new velocities that dissipate the maximum amount of energy without 
changing the convected distribution of mass) to the reader as they are similar to the 
BDF2 case. The generalization to higher order multi-step and diagonally implicit 
Runge-Kutta methods is also straightforward, although we have not tested them. 

2.3. Implementation details. 



2.1 and 2.4 



are easy 



The schemes of Definitions 
to implement once a solver for convex optimization is available. We experimented 
with our own implementation of a trust-region method described below, and with 
the convex optimization solvers cvxopt p] and Knitro 12 



In order to illustrate 



the key issues, let us consider the problem of minimizing the functional 
3 , 



F(z) 



N 

4=1 



m z. 



^"|2 



-f 



Af-1 

E 

i=l 



u 



H+1 



{Zt. 







(2.14) 



among all z G M such that z^+i > Zi for all i. The first term in (2.14) penalizes the 
difference between the positions Zi and the intermediate positions x", and since 1/r^ 
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is very large for small timesteps, it is natural to choose the intermediate positions as 
the starting vector for the minimization process. To have a bounded internal energy, 
we spread the x" slightly apart so that the minimal distance between neighboring 
positions of the starting vector is bigger than some number dmin. We define 

zf^ ■.^x\'^\{k + n), (2.15) 

where the negative displacements li are built up by 

for z = 2 to TV 

for j = i — 1 downto 1 

d = x^+i + Ij+i — (i" + Ij) — dmin 
if (i<0 

Ij — Ij + d (decrease Ij) 
else 

break (continue with outer loop) 

and a similar algorithm is used to obtain the r^. Note that this procedure preserves 
the ordering. In practice, only a few points i" need to be adjusted and most of the 



resulting positions z^ are equal to £". We use dmin := ^ mini(a:: 



2 ^^"^^H-^i+1 



X ■ 



From this initial guess, we proceed with a trust region Newton method to mini- 
mize the functional F{z). Our implementation is based on Chapter 4 of [14]. For 
a given z € M^, the gradient g = Vz;F(z) and Hessian H = W'^F{z) are trivial to 



compute for the first term in (2.141, and are easily assembled element by element 
(one interval [zi,Zi+i) at a time) for the second term. The Hessian is positive def- 
inite (as it is strictly diagonally dominant) and tridiagonal. At each step of the 
minimization algorithm, we seek the minimizer of the quadratic approximation 

Q(p):=F(z)+/p+ip^iJp « ^(z + p) 

over all p e M"* satisfying the trust region constraint ||Z3^^p|| ^ A. Here D is the 
diagonal scaling matrix with entries Da := ^ minimi — Zi_i, Zi+i — Zi} that prevents 
the positions Zi from crossing. The trust region radius A is not allowed to exceed 
one. By defining p := D~^■p^ we map the elliptical trust region to a sphere. We 
then consider the following modified problem: Minimize the functional 

Q(p):=F(z)+5^p+ip^iJp, (2.16) 

over all p e M^ with ||p|| < A. Here g := Dg and H := DHD. Notice that this 
minimization problem reduces to the previous one if I? = id. 

The main theorem governing the design of trust region methods is 

Theorem 2.8 (Trust region optimality criterion). The vector p^ e M.^ minimizes 
the quadratic functional (2.16) over all p G K^ with ||p|| ^ A i/ and only «/ p* is 



feasible and there exists a scalar A* ^ with the following properties: 

{H + X,id)p,^-g, (2.17) 

A,(A-||p,||) = 0, (2.18) 

{H + A* id) is positive semidefinite. (2.19) 



In our case, the matrix H is positive definite, so (2.19) is satisfied automatically. 



Moreover, equation (2.171 has a solution for any A* ^ 0, which eliminates the need 



to search for the most negative eigenvalue of H and also rules out the "hard case" 
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of More and Sorensen; see fl4|. We use a variant of the Cholesky factorization 
algorithm to find the Lagrange multipher A*. In pseudo-code it is given as 

Given A(o), A>0 
for fc = 0, 1, 2, . . . 

Factorize H + A^'^Hd = LL'^ 

Solve LL^pk = ~g and Lq^ — pk 

Define A(^+i) ^ AW _. f\M\f\M-^ 



..ml 

if A(*^) = and A^^+i) < 

return 
if |(|b,||-A)/A|< 10-12 

compute Pk+i 

return 
if A^'^+i) < 

set A^^+i) = 
end for 

This algorithm is equivalent to Newton's method for finding the zeros of the function 
A~i — ||p(A)||-i, which usually converges to machine precision in 4-6 iterations. 
We vary the size A of the trust region in the standard way J14l by comparing 
the actual reduction of the function to the predicted reduction of the quadratic 
model. For any given A, we use the Lagrange multiplier A.^ of the previous trust 
region step as starting value A*^"^ for the iteration above. Since we use the exact 
Hessian, the method converges quadratically (usually requiring 4-20 trust region 
steps). Moreover, as H is tridiagonal, the Cholesky factorization is of complexity 
0{N) only, and so the convex optimization procedure can be solved very efficiently. 



The minimization problem for the functional (2.9) is almost identical to the one 



we just described. The first term in (2.141 must be replaced by 

N 



E 



2^2 (^* XjAyXZj x^) ^^2 



[Zi Xj^)A.ij[Zj Xj) 2^^* •^iJ^Jil^i ■^j) 



The Hessian matrix for this part of the objective function is jpiA, which is again 
positive definite and tridiagonal, and so the problem can be solved efficiently. 



3. Numerical Experiments 
In this section, we report on numerical experiments we performed with the vari- 



ational particle methods VPSl and VPS2 introduced in Definitions 2.1 and 2.4 We 
studied both the porous medium and heat equations, as well as the isentropic and 
isothermal Euler equations, for various choices of parameters and initial data. 

3.1. Porous Medium and Heat Equations. For both the porous medium equa- 
tion (1.7) and the heat equation (|1.11[), solutions converge to self-similar profiles 



as i — >■ oo. These profiles can be computed explicitly, which makes it possible to 
measure the error of the numerical approximation exactly. 
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3.1.1. Porous Medium Equation. In the case of the porous medium equation the 
self-similar solution are called Barenblatt profiles and given by the formula 



where (s)-|_ := max{s,0}, and 



1/(7-1) 



(3.1) 



d 



d(7-l)+2' 



/3 = 



fc = 



/3(7^1) 
27 ■ 



see 
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The constant C is chosen in such a way that M{t) := Jg^ g^ (t, x) dx equals 
the total mass. Notice that M{t) is in fact independent of t > 0, and for simplicity 
we assume that M{t) ~ 1. Then gi*(i, •) converges to a Dirac measure as t — > 0. In 
the one-dimensional case under consideration, the constants simplify to 



1 



1' 



/3 = a, k 



7-1 
27(7 + 1)' 



and 




where r(z) := /^ i^ 



z ^ 0. To obtain an approximate solution of 



the porous medium equation (1.7) we used the variational particle schemes VPSl 
and VPS2 of the previous section with internal energy U{q) — Q^ l(j ~ 1). 



Figure IT] shows the solution of (1.71 with 7 = 5/3 at different times. We started 
out with asymmetric and discontinuous initial data 



q{x) 



0.5 if a; e (-1,0), 
0.25 if a: G (0,2), 
otherwise. 



We used the VPSl scheme with iV = 1000 points and r = 0.0016 to generate these 
plots. The running time was 8 seconds on a 2.4 GHZ laptop. Note that the solution 
at time T = 1.6 is close to a (shifted) Barenblatt profile. 

To estimate the accuracy of the VPSl scheme, we used initial data 



q(x) :-- 



if xe (-0.01,0.01), 
otherwise, 



(3.2) 



which approximates a Dirac mass, and computed the solution at time T = 10 for 



both the porous medium equation (1.7) with different values of 7 > 1, and for the 
heat equation, which corresponds to the limiting case 7 = 1. The particle mass was 
chosen as to = 0.001 and the timestep r = 0.01. The results are shown in Figure [2] 
Note that the solution is not continuously differentiable for 7^2: The contact 
angle is strictly positive for 7 = 2, is vertical for 7 > 2, and vanishes for 7 < 2. For 
the heat equation, the self-similar solution is given by a Gauss function 



g*{t,x) 



1 



V47ri 



exp 



\xr 
At 



(3.3) 
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0.5 



0.0 
0.5 



0.0 
0.5 



0.0 
0.5 



0.0 
0.5 



0.0 
0.5 



T = 0.0 




T = 0.1 




T = 0.2 




T = 0.4 



T = 0.8 



T= 1.6 




Figure 1. Porous medium equation with 7 = 5/3. 



The following table shows the ^°°-error of the approximate solution: 



Problem 


Jf°°-error 


heat eq. 


9.45e-5 


7 = 5/3 


1.53e-4 


7 = 3 


8.63e-4 


7 = 5 


2.62e-3 



This error is computed, with n := T/t, according to the formula 



.if °° -error 



max 



^i+1 



.(lO,i(xf 



^i+l. 



(3.4) 
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0.25 

0.00 
0.25 

0.00 
0.25 

0.00 
0.25 

0.00 



1 1 1 1 1 1 1 1 1 1 

heat 


1 
equation 


1 1 1 1 1 1 1 1 1 1 

7 


1 
= 5/3 






7 


= 3 




1 1 1 1 1 1 1 1 1 1 


= 5 

1 



-6-5-4-3-2-1012345 



Figure 2. Porous medium/heat equation at time T = 10.0. 



Since the Barenblatt profile has steep gradients for large exponents 7, the infinity 
norm of the error grows as we increase 7. On the other hand, the support of the 



solution of (1.7) decreases as we increase 7. For the heat equation, which has an 



infinite speed of propagation, the support is unbounded. The two effects largely 
cancel each other out when we estimate the ^-"^-error as follows: 



Problem 


^"'^-error 


heat eq. 


l.lle-3 


7-5/3 


1.08e-3 


7 = 3 


9.62e-4 


7-5 


7.26e-4 



Here the error is computed, with n :— T/r, according to the formula 



N-l 



^^-eiTOT := y. 



i=l 



• \^^^ 2 y^i ~^ ■^i+1. 



{xl_,-xf). (3.5) 



To compute the convergence rate of the VPSl method, we computed the solution 



at time T = 10 of the porous medium equation with initial data (3.2) and 7 = 5/3, 
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for several choices of particle mass m and timestep r. We have 







Error at 












m 


T 


x = 


Rate 


^°° -error 


Rate 


.if ^ -error 


Rate 


0.01 


0.1 


1.15e-3 




1.15e-3 




6.79e-3 




0.004 


0.04 


5.21e-4 


0.86 


5.21e-4 


0.86 


3.35e-3 


0.77 


0.001 


0.01 


1.53e-4 


0.88 


1.53e-5 


0.88 


1.08e-3 


0.82 


0.0004 


0.004 


6.73e-5 


0.90 


6.73e-5 


0.90 


4.98e-4 


0.84 


0.0001 


0.001 


1.90e-5 


0.91 


3.34e-5 


0.50 


1.47e-4 


0.88 



Thus, the VPSl method performs slightly worse than a first order method for these 
initial conditions. The last data point for the ^°°-norni is anomalous because the 
location of the largest error jumps from the center a; = to the boundary of the 
support of the solution, which is not a smooth transition. 

In this experiment, the initial data was rather singular. To estimate the effect of 
this irregularity, we also performed a computation starting off from the Barenblatt 
profile at time t = 1, which is much more regular, and continued to i = 2. More 
specifically, we partitioned the support [— C/vk, C/vfc] of g*{l,-) into intervals 
li = [xi-i,Xi] such that Jj g^{l,x)dx = 1/N for all 1 ^ i ^ iV, and define Xi to 
be the center of mass of g*(l, •) over /.;. Again we used 7 — 5/3. We have 







Error at 












m 


T 


x^O 


Rate 


.if°°-error 


Rate 


^^-error 


Rate 


0.01 


0.05 


4.32e-4 




7.99e-4 




1.58e-3 




0.004 


0.02 


1.88e-4 


0.91 


5.69e-4 


0.37 


7.88e-4 


0.76 


0.001 


0.005 


5.00e-5 


0.96 


2.89e-4 


0.49 


2.38e-4 


0.86 


0.0004 


0.002 


2.04e-5 


0.98 


1.76e-4 


0.54 


1.03e-4 


0.91 


0.0001 


0.0005 


5.17e-6 


0.99 


8.03e-5 


0.57 


2.80e-5 


0.94 



In the interior of the support of the solution, the method converges at first order, 
but near the boundary the order deteriorates. The .if ^-error involves a mixture of 
the two convergence rates, and is therefore slightly worse than first order. 

We also tested the convergence rates of the VPS2 method. For initial data of 
the form (3.2) with 7 := 5/3, we chose TV -j- 1 knots Xi uniformly spaced over the 
interval [—0.01,0.01] and set rrii = 1/N for all 1 ^ i ^ iV. The errors are smaller 
than they were for VPSl, but the method still converges at first order only due to 
the discontinuity in the initial density. We have 
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To estimate the influence of the singularity of the initial data, we tried the VPS2 
method starting again off from the Barenblatt profile at time t = 1 and continuing to 
t — 2. We partitioned the support [—C/Vk, C/Vk] =: [a, b] into equal subintervals 
with knots Xi = a + i{b — a)/N for all 1 ^ z ^ A^ and set 



TO,; 



(l,x)dx for ah 1 s$ i ^ A^. 



(3.6) 
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We used 7 = 5/3. When comparing the result to g*{2, •), we have 
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As expected, we obtain second order convergence throughout most of the interval 
when we start with continuously differentiable initial data. A few points next to 
the boundary of the support, however, seem to converge at a lower rate of 3/2. 
This raises the question of what happens if the initial data is continuous but not 
'^^^. The Barenblatt profile with 7 = 3 has a square-root singularity at the edges 
of its support. Repeating the above procedure for this case, we obtain 
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It appears that the infinite slope of the exact solution near the boundary of the 
support creates large errors there that dominate the overall .jSf^-error of the method. 
The situation can be improved if we redistribute the mass to better resolve the 
solution near these endpoints. Specifically, we choose the knot positions as 

^ n) ^ 

with weight function 



/ 



for all s^ i s^ A^, 



/(^) 



yjl - y2 ^y 



Jo 



for all X e [-1, 1], 



and we assign to each interval [xi^i,Xi) the mass rrii defined in (3.6). We have 
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In spite of the singularity in slope near the endpoints of the support, we achieve 
second order accuracy over enough of the interval to converge at second order in the 
.jSf ^-norm. The ^°°-norm is larger because the first and last intervals are smaller 
than the interior intervals, so their midpoints are closer to the singularity. The 
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first function we tried for the mass redistribution, namely f{x) 
difficulties for the convex optimization solver for N — 10000. 



Ix^, caused 



3.1.2. Heat Equation. Both schemes perform similarly for the heat equation ( |1.11[ ) 
as they did for the porous medium equation. Here we report only the VPS2 results. 
Using the internal energy U{g) — g\ogg and initial da ta p.2[ ), and comparing the 
solution at i = 10 to the Gaussian q^{10, •) defined in (3.3), we find 
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As before, we obtain first order convergence for the VPS2 scheme for discontinuous 
initial data. Note that we could obtain a more precise error estimate by comparing 
the numerical solution to the exact solution of ( |1.11[ ) with data (3.2 1, which can be 
computed explicitly. But the ^""-distance between this exact solution and g*(10, •) 
is around 7.4 x 10~^, which is small in comparison to the numerical errors. 

Finally, we use the VPS 2 scheme to evolve the heat kernel e*(l, •) to time t — 2 
via the heat equation ( |1.11[ ). As the support of gi*(l, •) is infinite, instead of fixing 
the knot positions first, we start by distributing the mass according to 



/ 



N 



N + 1 



for a\ll€i€N, 



E/ 



iV + l 
1=1 

with weight function 

f{x) :— q{x)q{l — x) and q{x) :— lOa;^ -|- x/lO. 
We then initialize the knot positions a;? for all 1 ^ i ^ A^ — 1 as 



-2erfc~\2s,) with 



s, := 



E^ 



This construction implies that rrii = J ^ g^,{l,x) dx for all 2 < i < A — 1. Since 
the mass of the first and last interval are small, any reasonable choice of xq and xn 



works well. We used 



3a;; 



^Xn anQ X j\j 



•jXj^_ 



2x%_2- We have 
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Thus, for the heat equation, even the ^°°-norm converges at second order as there 
is no singularity in the slope of the exact solution over the course of this simulation. 
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The ^^-norm converges slightly slower than second order due to the fact that as 
we add points, the support of the numerical solution grows. 

In conclusion, we have found our variational particle scheme to be an effective 
method of solving the porous medium and heat equations. Moreover, the second 
order version of the method does indeed give second order accuracy as long as the 
initial conditions are not too singular. 



3.2. Isentropic Euler Equations. As few smooth solutions of the isentropic Eu- 



ler equations (2.1) are known explicitly, we test our numerical methods for self- 
consistency using smooth initial data for a short enough time that shocks do not 
occur. To study their ability to capture shocks and rarefaction waves, we test our 
schemes using piecewise constant initial data, for which the Riemann problems 
can be solved analytically. In general, we find that shocks cause the second order 
method to perform as a first order method, and the first order method to perform 
slightly worse. We mention in passing that these schemes also work well for the 
isothermal equations, but we did not include examples for the sake of brevity. 



3.2.1. Smooth solutions. We begin by testing the convergence rates of the VPSl 
and VPS2 schemes for the isentropic Euler equations with initial data 

3/ a;2\ 

q{x) = - I 1 I and ii{x) — for all a; e M, 

8 V 4 /^ 

for 7 = 5 and 7 = 7/5. The solutions are shown in Figure Is] at times i = 0,l,2,3,4. 
For the VPSl scheme, we choose initial positions a;^ such that 

r^° i -1/2 

/ 5(x) dx = '— for aU 1 < i < A^ . 

For the VPS2 scheme with 7 = 5, we choose initial position x? uniformly spaced in 
the support [—2,2] of the initial density, and define masses 



g{x, 0) dx for all 1 s^ i sC iV. (3.7) 



For the VPS2 scheme with 7 — 7/5, it is necessary to distribute more points near 
the ends of the support to achieve full second order accuracy. We found that 

X, :== 2/(-l + 2i/N) for aU < i < A, 

with weight function 

f{x) := ^ for all x E [-1, 1], 

v/l - j/2 dy 
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evolution of density 
t=0 



evolution of velocity 




Figure 3. Solution of isentropic Euler equations with a parabolic 
initial distribution of mass, starting from rest. {N — 800, VPS2) 



and with rrii defined in (3.7) works well. We used the VPS2 solution with N — 4000 



and T — 0.004 as the "exact" solution to compute the errors for both schemes. 
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The errors reported for the VPSl scheme are estimates of the ^°°-error of q alone, 
measured at midpoints of the numerical solution against the "exact" solution, in- 
terpolated linearly between its own midpoints. The errors reported for the VPS2 
scheme involve the Wasserstein distance of the densities and the difference in ve- 
locities. Notice that the finiteness of the total energies only implies that 



(^^'^{^,0^) and 



e.if2(M'',g'^^'=*), 
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SO in general the two velocities are in different spaces. Let r" and r'^^^'^^ denote the 
inverse distribution functions of the measures gi"£^ and g™a'Ct£i j.ggp ^g explained 
we then have r"#(l[o^i]£^) = (f^S,^, which implies that 



in Section 



2.2 



That is, the pull-back uo r" of the velocity field m" under the map r" is an element 
in ^^([0, 1]). The same is true for u^^^^* o r™^'^*. We therefore define 

V ^ "'[0,1] 

where i" = 4 is the final time. Since convergence in the Wasserstein space implies 
only weak* convergence in the sense of mesaures (see [TIITt]), our error estimate is 
weaker that one that involves the .jSf^-norm. The ^°°-error of the densities of the 
VPS2 solutions are about three times larger than the i<^vv-errors when 7 = 5, and 
about three times smaller than the E'w-errors when 7 = 7/5. 

In the following three sections, we study the performance of our scheme on 
problems involving various combinations of shocks and rarefaction waves emanating 
from an initial discontinuity at the origin. In the presence of shocks, it no longer 
makes sense to use the ^°°-norm to measure the error. 

3.2.2. Shock/Shock. We use 7 ~ 5/3 and initial data 

[(0.25,1) if.Te(-2,0), 
(g,w)= < (0.25,0) ifxe(0,2), (3.9) 

[ (0, 0) otherwise. 

Note that this amounts to solving three Riemann problems: one at x = 0, and two 
others at a; = ±2. Figure |4] shows the evolution at different times, computed via a 
first order version of the VPS2 method (with variable mass 

"^^ ^= Hl(^) (' " ^) '" foralll^^^iV, 

but a first order time-stepper), with A^ = 1000 and timestep r — 0.01. We will refer 
to this hybrid method as VPSla. Since there cannot be any shocks connected to the 
vacuum, the solution immediately forms two rarefaction waves at the boundaries of 
the support that grow in width in a self-similar fashion. The Riemann problem at 
a; = evolves into a self-similar, double shock structure with constant intermediate 
state {gm,Um) — (1.16641,0.5) with large density. For simplicity, we only plot the 
density here, omitting the velocity. As time moves on, the different waves eventually 
interact and the shock strengths decrease. At time t — 7, the solution has developed 
into a continuous profile (not unlike the Barenblatt profiles of the porous medium 
equation), which seems to evolve in a self-similar way while moving to the right 
with a constant background speed. 

The VPSl and VPS2 schemes give similar results to this hybrid method, but the 
former is less accurate inside the rarefaction waves (due to lack of resolution) while 
the latter suffers from small oscillations near shocks. These oscillations are rather 
interesting as they do not grow exponentially in time nor prevent the scheme from 
converging as we refine the mesh; see Figure [5J Instead, they arise because once we 
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Figure 4. Isentropic Euler equations with 7 = 5/3. 



choose a spatial discretization {N together with the mass rrii of each interval) , the 
time-stepper is effectively solving a coupled system of masses connected by springs. 
Indeed, except for the step in which mass is re-distributed when the particles cross, 
we are solving Lagrange's equations of motion for the particle trajectories Xi(t), 
with the kinetic and potential energy of the system defined in terms of the Xi and 
Ui = Xi in the obvious way (depending on the scheme). If we were to take the 
time-step r to zero holding N constant, the particles would never cross during the 
transport step and the scheme would convert the energy that is supposed to be 
dissipated by the shock into lattice vibrations. However, by simultaneously refining 
the mesh as we decrease the time-step (with r ^ N~^), the backward-Euler and 
BDF2 methods seem to damp out these vibrations without significantly smoothing 
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Figure 5. Comparison of the VPS2 and VPSIa results with the 
exact solution of the Riemann problems (dashed lines). 



the shocks, which is somewhat amazing. We emphasize that although the VPSIa 
scheme does not suffer from oscillations for the timestep used in Figure [5] they will 
eventually appear if we take r — >■ holding N constant. 

Next we wish to use the exact solution of the Riemann problems to compute the 
errors of the schemes and check convergence rates. For sufficiently small t > 0, the 
exact density is given by 



£i(a;,i) 



e + 1 



i/s 



Qi 

Qm 
Qr 



{Ur + Q^j) -{x- Xr)/t 



1/0 



if x e (a;i,X2), 

if a; e (a;2,a;3), 
if a; e {xz.Xi), 
if a; e (0:4, 2:5), 

if a; e {x^.xq), 

otherwise, 
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Figure 6. Energy (solid: total, dashed: kinetic, dotted: internal). 



with constants 

xi=xi+t{ui- g1), X3^tsi, 

see pi. For the case under consideration, we have 
XI = -2, gi = 0.25, 

Xr =2, Qr — 0.25, 

which gives an intermediate density p„j — 1.16641, and shock speeds si = 0.36360 
and Sr = 0.63640. The exact velocity profile can be computed as well, but we refer 
the reader to the literature. Using a final time of T = 1.6, we obtain 



X5 = Xr +t{Ur - 9gr), 




xe ^ Xr + t{ur + g^)'^ 




ui = 1, 


(3.10) 


Ur =0, 


(3.11) 
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Here W stands for the Wasserstein distance between the numerical and the exact 
densities, while ^^ indicates their ^^ -difference. The functional Ew was defined 
in (3.8), and we write i^tot 



:= E" 



TT'cxact I 

" Etot I 



for the difference in the total energies. 
The .if °°-error of the densities does not approach zero because of the shocks. We see 
that in spite of the high frequency oscillations, the VPS2 scheme is more accurate 
than the VPSla method, with the density converging at first order and the velocity 
converging somewhat slower. For both schemes, the total energy E^^^. at the final 
time converges to fi't™'^' from above as A^ — > oo and r — >■ 0, thus the exact solution 
dissipates energy (slightly) faster than any of the numerical solutions. 

In Figure l6J we plot the energy of the numerical solution (beyond the point that 
we are able to compute the exact solution). As long as the solution is discontinuous. 



VARIATIONAL PARTICLE SCHEME 



29 



the total energy decreases. Moreover, the total energy is a linear function at the 
beginning of the evolution (up to time t — 2, say,) since the strength of the two 
shocks remains constant. During this time, the internal energy increases along with 
the width of the intermediate state (pm, Um), which carries more energy because of 
its high density. As the shocks interact with the rarefaction waves at the boundary 
of the support, their strengths decrease and finally vanish (at about t = 6.7), after 
which the total energy remains essentially constant. This is in agreement with the 
theory since for smooth solutions of the isentropic Euler equations, the total energy 
is a conserved quantity. Our schemes capture this behavior remarkably well. 



3.2.3. Shock/Rarefaction. We also considered the case 7 — 5/3 with initial data 

f (0.5,0) ifxe(-l,0), 
(^,m)= < (0.25,0) ifxe(0,2), (3.12) 



(0.5,0) ifxe(-l,0), 
(0.25,0) ifxe(0,2), 
(0, 0) otherwise. 



The exact solution of the Riemann problem at a; = is given by a pattern involving 
a shock followed a rarefaction wave. Figure [7| shows the approximation at different 
times. The computation was done using VPSl with m = 0.001 and r — 0.01. 

Figure ^ shows the approximate solution of ( |2.1| for the same parameters and 
initial data p.l2 |, at time T — 0.6 for different values of 7. For sufficiently small 
i > 0, the exact density is given by the following formula: 



g{t,x) 



i-ui + ef ) + {x- Xl)/t 



i/e 



Qi 



[ui + gf ) - x/t 

e + 1 



i/e 



[Ur + g^) ^ [x - Xr)/t 



1/e 



if x e (a;i,a;2), 
if a; e {x2,xs), 

if a; e (2:3,2:4), 

if a; e (a;4,a;5), 
if a; e (a;5,a;6), 

if a; e (a;6,a;7), 
otherwise. 



(3.13) 



with constants 

a:i = xi +t{ui - g1), 
X2 = xi +t{ui +0g^), 



X3 ^t(ui -Oqi), 

Xi = t{Uyn ~ OqI^), 
X^ — ^^r. 



Xq 
XT 






t{Ur + qI), 



see pj. For the case under consideration, we have 



XI = -1, 


QI = 0.5, 


ui =0, 


Xr = 2, 


Qr = 0.25, 


Ur = 0, 



which implies an intermediate density Qm — 0.3601 and velocity Um — 0.0823, and 
a shock speed s^ = 0.3636. The exact velocity profile can be computed as well, but 
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Figure 7. Isentropic Euler equations with 7 = 5/3. 



we refer the reader to the hterature. The fohowing table shows the ^^-error: 



Problem 


^^-error 


7 = 5/3 


4.46C-3 


7 = 3 


5.57e-3 


7 = 5 


3.87e-3 



Notice that the rarefaction waves in (3.13) that connect the profile to the vacuum, 



are convex for 7 < 3, linear if 7 = 3, and concave if 7 > 3. To get a convergence 



rate, we computed the approximate solution at time t = 0.6 with initial data (3.12) 
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Figure 8. Isentropic Euler equations at time T — 0.6. 



and 7 — 5/3 for different values of particle mass/timestep. We have 



m 


T 


^^-Error 


Rate 


0.01 


0.1 


2.146C-2 




0.005 


0.05 


1.397C-2 


0.620 


0.001 


0.01 


4.464e-3 


0.709 


0.0005 


0.005 


2.615e-3 


0.772 


0.0001 


0.001 


7.318C-4 


0.791 



As with VPSla in the previous section, the rate is somewhat less than one. We did 
not study the performance of the VPSla or VPS2 schemes on this initial data or 
implement the other measures of error in our VPSl code. 



3.2.4. Rarejaction/ Rarefaction. Finally, we use the VPSl method to compute the 
approximate solution of (2.1 ) for 7 = 5/3 and initial data 



{q^u) 



(0.25,-0.5) 

(0.25,0.5) 

(0,0) 



ifxe (-2,0), 
ifxe (0,2), 
otherwise. 



Figure |9] shows the result at different times. The fluid splits into two parts that 
travel in opposite directions. Between the blocks, two rarefaction waves form that 
are separated by vacuum. In our approximation, the density is defined as particle 
mass divided by the specific volume, which is the distance between neighboring 
particles. We therefore do not get a perfect vacuum. As the two blocks move further 
apart, however, the density between them does approach zero. By construction, 
the density can never become negative. 
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Figure 9. Isentropic Euler equations with 7 = 5/3. 
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Appendix A. Porous Medium Equation 
We derive the Euler-Lagrange equations for the minimization problem 

£)"+! := argmin \ ^W{g'\ g)'^ + U[g] : g € ^{W^) \ 



(A.l) 
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to show that it indeed gives an approximation for the porous medium equation. 
More precisely, we will show that (A.ll can be interpreted as a backward Euler 
method applied to Darcy's law u — —'VU'{g) in a Lagrangian formulation of the 
problem. We recall first that if g" and g are Lebesgue measurable functions, then 
there exists a Borel map r : M'* — > M'^ defined £»"-a.e., with the property that 

W(e",g)2=/" |r(a;)-xpp"(x)dx and r#(g"£'*) = £.£'*. (A.2) 

The map r is called an optimal transport map. It is the gradient of a lower semi- 
continuous, convex function and invertiblc g"-a.e. The second identity in (A.2 ) now 
implies that for all test functions (p E '^h{^'^) we have 

ip(r{x))g^{x)dx ~ / ip{z)g{z)dz — / ip{r{x))g(^r{x)) dct Dr{x) dx. 

Here we used the change of variables formula, which can be justified because r is a 
monotone function. Since r is invertible g"-a.e. we conclude that 

g"{x) = g{rix)) det Dr{x) for a.e. x G E'', 

which allows us to express the internal energy U[g\ in the form 

g"{x) 



U[g] 



U{g{z))dz^ / U 



det_Dr(a;) dx. 



(A.3) 



^detDv{x) 

That is, we can consider the internal energy as a functional of r instead of g. 

Let now £i"+^ be the minimizer of (1.12) and let r"+^ denote the optimal trans- 
port map that pushes £»"£'' forward to the measure g'^'^^2,'^. For any smooth vector 
field C : W^ — > W^ and any e > we now define the functions 

Y, := (id + eC) ° 1-"+^ and g^^"^ := r^#(e"£''). 

Note that the map r^ is not an optimal transport map if e > 0, but we have 

W(ei",e^)2 < /" |re-idpg"dx foraUe^O. 

This implies the estimate 

limsup^(^"'^-)^-^(^"'^"")%2/ (Cor"-).(r"--id),".. (A.4) 

£^0 e jRd 



C • ((r"+i - id) o (r"+i)-i)e"+i dz. 



On the other hand, a straightforward computation using (1.2 ) and (A.3 1 shows that 



e^O e 



U' 



g" 



U 



det Dr"+i 



X tr 



detL»r"+i/ deti^r'^+i 
({DC)or'^+A detDr''+^dx 

Jr'^ 
We used the change of variables formula again. Now (A.l) implies that 

Os=:- / C- ((r"+^-id)o(r"+i)-i)(?"+idz- / P{g''+^)V ■ (dz 



(A.5) 
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for all test functions C,. Note that by changing the sign of C we obtain the converse 
inequality, so we actually have equality here. Since C was arbitrary we obtain 



-id 



(r"+i)- 



„n+l 



VP(e"+^) = 



in the sense of distributions. If we define u""""^ := T~^(r"+^ — id)o(r""'"^)~-'^, then we 
obtain Darcy's law u""*"^ — —WU'{g"~^^) a.e. in {g"^^ > 0}. Since the new density 
p"+^ can be computed according to the formula p" = io"'^^ or"+^) det£'r"+^ a.e., 



we conclude that the minimization problem (1.121 in fact amounts to a backward 



Euler method for the transport map r"+^ pushing p'^fi"^ forward to g^~^^£'^ 



Appendix B. Isentropic Euler Equations 

The time discretization of Definition [Ll] is motivated by Dafermos' entropy rate 
criterion [6]: dissipate the total energy as fast as possible. This objective is achieved 
through two different mechanisms. First, the kinetic energy is minimized by replac- 
ing the given velocity by the corresponding optimal transport velocity. Second, the 
internal energy is minimized by solving an optimization problem similar to the one 
in (1.12) for the porous medium equation, but with a different penalty term. The 



intuition for the latter is the following: We think of the fluid as a collection of parti- 
cles, each one determined by a position and a velocity. The particles prefer to stay 
on their free flow trajectories (minimizing the acceleration), but the pressure from 
the surrounding particles forces them to deviate from their characteristic paths. 

To measure the acceleration we introduce a functional similar to the Wasserstein 
distance. Consider first a single particle that initially is located at position xi G K'^ 
with initial velocity ^i € M''. Specify some timestep r > 0. During the time interval 
of length T we allow the particle to move to a new position X2 G K^ and to change 
its velocity to ^2 G K''. What is the minimal acceleration required to achieve this 



change? We are looking for a curve c: [0,t] 



such that 



(c,c)(0) = (a;i,^i) and (c, c)(t) = (0:2,6). 

that minimizes the .if^([0, r])-norm of the second derivative c along the curve. 
This minimization problem has a unique solution; the minimizer is given by a cubic 
polynomial; and the minimal (averaged) acceleration is given by the formula 



\'c{s)\^ds^ 12 



X2 - Xi 



6 + 6 



6-6 



Consider now a fluid characterized by a density g and a velocity field u. Assume 
that we are given a transport map r: M'* — > M'' that is invertible g-a.e. and deter- 
mines where each particle travels to during a time interval of length r > 0. Which 
velocity field u minimizes the total acceleration cost, defined as the integral 



<,, 



id 



-|uor- 



gdx? 



(B.l) 



Note that since r is assumed to be invertible, there is only a single particle located 
at position r(a;) for g-a..e. x G K'', and so the velocity of the transported particles 
can be described by an Eulerian velocity field u. Let u+ denote the minimizer of 
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the functional s^r.v (we refer the reader to [8j for further details). For any smooth 
vector field C: M'' — > W^ and e > define u^ := u+ + eC^. Then we compute 



lim — ■ — 

£^0 



12 



id 



2u^ 



(C o r)gdx, 



which vanishes because u"*" is a minimizer. Since C, was arbitrary, we obtain that 

. 3 



2r 



r — (id + ru) 



(B.2) 



g-a,.e. Using this identity in (|B.1|) we have 



3 

4^2 



/ |r-(id 



Tu)\'^Qdx =: 



Assume now that the push-forward of gH'^ under the map id + ru is absolutely con- 
tinuous with respect to the Lebesgue measure and define g£'^ := (id -f Tu)#(p£'^). 
If we fix a density g^ £ ^{M.'''), then there are many transport maps r : M"^ — > M'^ 
that push g£'^ forward to the measure g'^2,'^. Minimizing in r, we obtain 

inf {#;[r]2: r#(e£'') = g+Q'^} - ^W{g,g+r, 

so the minimal total acceleration cost is proportional to the Wasserstein distance 
between g'^ and the density g obtained by freely transporting the initial density g 
along the initial velocity field u. We refer again to M for further details. 

Note that Step 2 of the time discretization in Definition |1.1| amounts to minimiz- 
ing the internal energy U[g\ while penalizing the deviation of particles from their 
characteristic trajectories (which requires accelerating the particles). Computing 
the Euler-Lagrange equations for the minimization problem (1.12), as we did above 



for the time step for the porous medium equation, we find that the map 

2r2 



i.n+1 



= id 



3 



-VU'ig 



n+l^i 



is an optimal transport map between g"+^£'* and ^"£'*, thus invertible gi"+^-a.e. 



The update of the velocity (1.18) follows from formula (B.2) if we set 

2r2 



u := u", u+ .— L 

Indeed we have that 

3 / /. . 2r2 



,n+l 



id 



-WU'ig 



?i+i> 




-id 



WU'ig 



n+l 




(id 



Vf/'(e"+i) 



rVC/'(e"+i) 



£i"+^-a.e. The only caveat is that (B.2) was derived under the assumption that the 



map r is essentially invertible. But this condition is satisfied because u" was chosen 
in such a way that the map id -I- ru" is an optimal transport map pushing p"£'^ 
forward to g"'£'' := (id-|-ru")#(£i"£'^); see Step 1 in Definition |l.l[ As mentioned 
there, this step decreases the kinetic energy in an optimal way. 
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